In the implicit method called alpha schema, used to solve non-stationary problems, our differential operator describing the "physics" of the simulated phenomenon is applied to the linear combination of the previous and next state values:
\( \frac{u_{t+1}-u_t}{dt} - \mathcal{L}(\alpha u_t+(1-\alpha)u_{t+1})) = {f_{t+1}} \).
If the differential operator is linear, then we can break it into two terms
\( u_{t+1} -dt \alpha \mathcal{L}(u_{t+1}) = u_t+dt{f_{t+1}}+ dt(1-\alpha)\mathcal{L}(u_t) \)
leaving our "unknown" on the left-hand side \( u_{t+1} \).
Note that for the value of alpha = 1, this scheme is an implicit Euler scheme (the so-called backward Euler), for the value of alpha = 0 it is an explicit scheme (the so-called Euler scheme) and for the value of alpha = 1/2 it is the Crank-Nicolson scheme (which is also implicit).
The Euler scheme was first described in 1768 in a textbook by Leonhard Euler. The Crank-Nicolson diagram was invented in 1947 by Phyllis Nicolson, a female British mathematician, and John Crank, an American mathematician and physicist [1], [2].
We use the isogeometric finite element method.
To develop a solver that uses the method implicitly in the isogeometric finite element method, we need to transform a strong formulation into a weak formulation.So we multiply our equation by the test functions \( (u_{t+1},w) -dt \alpha(\mathcal{L}(u_{t+1}),w) = (u_t,w)+dt({f_{t+1}},w)+ dt(1-\alpha)(\mathcal{L}(u_t),w) \).
We will use a linear combination of the B-spline function to approximate the state of our system at a given moment in time. For this purpose, we select the basis of two-dimensional B-spline functions, defining the node vectors in the direction of the coordinate system axis. To establish attention, we can choose a two-dimensional basis of the second order B-spline function
\( \{B^x_{i,2}(x)B^y_{j,2}(y)\}_{i=1,...,N_x;j=1,...,N_y} u_{i,j } \).
They will be used to approximate the simulated scalar field of the current time instant \( u(x,y;t+1)\approx \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } B^x_i(x)B^y_j(y) \)
Similarly, for testing, we will use the B-spline base functions:
\( w(x,y) \in \{ B^x_k(x)B^y_l(y) \}_{k=1,...,N_x;l=1,...,N_y} w^{k,l } \)
Assuming that the differential operator describing "physics" is linear, our equation now looks like this: \( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } (B^x_i(x)B^y_j(y)-dt \alpha \mathcal{L}(B^x_i(x)B^y_j(y)),B^x_k(x)B^y_l(y))=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \). For example, for the problem of heat transport we have
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ (B^x_i(x)B^y_j(y)-dt \alpha \left(\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial x^2}+\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial y^2 }\right),B^x_k(x)B^y_l(y))=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \)
Thanks to the weak formulation, we can integrate by parts \( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ \left( (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) - dt \alpha(\frac{\partial (B^x_i(x)B^y_j(y))}{\partial x},\frac{\partial (B^x_k(x)B^y_l(y))}{\partial x}) -dt* (\frac{\partial (B^x_i(x)B^y_j(y))}{\partial y},\frac{\partial (B^x_k(x)B^y_l(y))}{\partial y }) \right)=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \)
Due to the structure of the tensor product of the B-spline function and due to the fact that the derivative in the y direction of the B-spline in the x-axis direction gives 0 (because these functions are constant in the y-axis direction) and similarly for the derivative in the y-direction of the B-spline in the x-axis direction we have
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \\ \left( (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) - dt \alpha (\frac{\partial B^x_i(x)}{\partial x}B^y_j(y),\frac{\partial B^x_k(x)}{\partial x}B^y_l(y)) -dt* (B^x_i(x)\frac{\partial B^y_j(y)}{\partial y},B^x_k(x)\frac{\partial B^y_l(y)}{\partial y }) \right)=RHS \\ \forall k=1,...,N_x; l=1,...,N_y \)
Note that our system of equations to solve is \( \left(M_x \otimes M_y-dt \alpha S_x \otimes M_y -dt* M_x \otimes S_y\right) u^{t+1} = F(t) \)
where
\( \{ M_x \otimes M_y \}_{i,j,k,l}= \int B^x_i(x)B^x_k(x)B^y_j(y)B^y_l(y) = \int B^x_i(x)B^y_j(y)B^x_k(x)B^y_l(y) = {\bf M}_{i,j,k,l } \) is a mass matrix which is a Kronecker product of two one-dimensional mass matrices,
\( \{ S_x \otimes M_y\}_{i,j,k,l} = \int \frac{\partial B^x_i(x)}{\partial x}\frac{\partial B^x_k(x)}{\partial x}B^y_j(y)B^y_l(y) = \int \frac{\partial B^x_i(x)}{\partial x}B^y_j(y)\frac{\partial B^x_k(x)}{\partial x }B^y_l(y) \) is the Kronecker product of the one-dimensional stiffness matrix and the one-dimensional mass matrix, and
\( \{ M_x \otimes S_y \}_{i,j,k,l} = \int B^x_i(x) B^x_k(x) \frac{\partial B^y_j(y)}{\partial y} \frac{\partial B^y_l(y)}{\partial y } = \int B^x_i(x)\frac{\partial B^y_j(y)}{\partial y}B^x_k(x)\frac{\partial B^y_l(y)}{\partial y } \) is the Kronecker product of the one-dimensional mass matrix and the one-dimensional stiffness matrix.
Each of these matrices can be factored in linear time using the variable-directional solver algorithm. However, it is no longer possible to factorize their sum in linear time. This is possible only when we introduce a time step scheme that allows the matrix to be split into sub-matrices in time sub-steps so that the factorization cost remains linear.
The Peaceman-Reachford scheme allows you to split the time step into two sub-steps \( \left(M_x \otimes M_y-dt \alpha S_x \otimes M_y\right) u^{t+1/2} = F(t+1/2)+\left(dt (1-\alpha) M_x \otimes S_y\right) u^{t} \)
\( \left(M_x \otimes M_y-dt \alpha M_x \otimes S_y\right) u^{t+1} = F(t+1/2)+\left(dt (1-\alpha) S_x \otimes M_y\right) u^{t+1/2 } \)
where we can merge the left-hand side matrices into a single matrix with a Kronecker product structure that can be factored in linear time using the variable direction solver algorithm
\( \left(M_x-dt \alpha S_x \right)\otimes M_y u^{t+1/2} = F(t+1/2)+\left(dt (1-\alpha) M_x \otimes S_y\right) u^{t} \)
\( M_x \otimes \left( M_y-dt \alpha S_y\right) u^{t+1} = F(t+1/2)+\left(dt (1-\alpha) S_x \otimes M_y\right) u^{t+1/2 } \).
Finally, let us take a look at what the right-hand side operator looks like.
The right-hand side function here represents the changes caused by the force acting on the system during the time step, plus the changes caused by the physics of the phenomenon in the previous time step \( (u_{t+1},w) -dt (1-\alpha)(\mathcal{L}(u_{t+1}),w) = (u_t,B^x_k(x)B^y_l(y))+dt({f_{t+1}},B^x_k(x)B^y_l(y))+ dt(1-\alpha)(\mathcal{L}(u_t),B^x_k(x)B^y_l(y)) \)
In particular, our right side is not a bitmap, but a projection of the sum of three elements:

1. The state of our system in previous time step \( (u_t,w)= \sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y)) \) (also multiplied by the test functions and integrated over the domain). Of course, we also use a linear combination of B-spline basis functions to represent the state of our system in the previous time step \( u(x,y;t)=u_t=\sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } B^x_i(x)B^y_j(y) \).
2. Changes induced by the "physics" of the simulated phenomenon during the time step. These changes are calculated by applying the differential operator representing the modeled phenomenon to the state of the system at the previous moment. The state of the system is of course represented by a linear combination of the B-spline function. So the differential operator is applied to the B-spline function \( (dt (1-\alpha) \mathcal{L}(u_t),w)=\sum_{i=1,...,N_x;j=1,...,N_y} dt*u^{t}_{i,j } ({\cal L} (B^x_i(x)B^y_j(y)),B^x_k(x)B^y_l(y)) \). For example, for the problem of heat transport we have \( \sum_{i=1,...,N_x;j=1,...,N_y} dt (1-\alpha) u^{t}_{i,j } (\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial x^2}+\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial y^2},B^x_k(x)B^y_l(y)) \). Thanks to the weak formulation, we can integrate by parts \( \sum_{i=1,...,N_x;j=1,...,N_y} dt (1-\alpha) u^{t}_{i,j } \left( (\frac{\partial (B^x_i(x)B^y_j(y))}{\partial x}B^y_j(y),\frac{\partial (B^x_k(x)B^y_l(y))}{\partial x}B^y_l(y))+(B^x_i(x)\frac{\partial (B^x_i(x)B^y_j(y))}{\partial y},B^x_k(x)\frac{\partial (B^x_k(x)B^y_l(y))}{\partial y})\right) \) which due to the structure of the tensor product of the B-spline function, and due to the fact that the derivative in the y direction of the B-spline function of the variable x gives 0, and the derivative in the x direction of the B-spline function of the variable y gives, we get \( \sum_{i=1,...,N_x;j=1,...,N_y} dt (1-\alpha) u^{t}_{i,j } \left( (\frac{\partial B^x_i(x)}{\partial x}B^y_j(y),\frac{\partial B^x_k(x)}{\partial x}B^y_l(y))+(B^x_i(x)\frac{\partial B^y_j(y)}{\partial y},B^x_k(x)\frac{\partial B^y_l(y)}{\partial y})\right) \). From the "matrix" point of view, this term gives us \( dt (1-\alpha) (S_x \otimes M_x +M_x\otimes S_y) u_t \).
3. Changes caused by a force acting on the system during the time step \( (f_{t+1},w)= (f_{t+1}(x,y),B^x_k(x)B^y_l(y)) \).